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1 Introduction 

Passive remote sensing of the atmosphere is used to determine the atmospheric state. A radiometer measures 
microwave emissions from earth's atmosphere and surface. The radiance measured by the radiometer is 
proportional to the brightness temperature. This brightness temperature can be used to estimate atmospheric 
parameters such as temperature and water vapor content. These quantities are of primary importance for 
different applications in meteorology, oceanography, and geophysical sciences. Depending on the range in the 
electromagnetic spectrum being measured by the radiometer and the atmospheric quantities to be estimated, 
the retrieval or inverse problem of determining atmospheric parameters from brightness temperature might be 
linear or nonlinear. In most applications, the retrieval problem requires the inversion of a Fredholm integral 
equation of the first kind making this an ill-posed problem. The numerical solution of the retrieval problem 
requi res the transformation of the continuous problem into a discrete problem. The ill-posedness of the 
continuous problem translates into ill-conditioning or ill-posedness of the discrete problem. Regularization 
methods are used to convert the ill-posed problem into a well-posed one. 

I n this paper, we present some results of our work in applying different regularization techniques to 
atmospheric temperature retrievals using brightness temperatures measured with the SSM/T- 1 sensor. Sim- 
ulation results are presented which show the potential of these techniques to improve temperature retrievals. 
In particular, no statistical assumptions are needed and the algorithms were capable of correctly estimating 
the temperature profi le corner at the tropopause independent of the initial guess. 


2 Radiative Transfer Theory in the Microwave Region 

Radiative transfer theory describes the intensity of radiation propagating in a general class of media that 
absorbs, emit, and scatter the radiation [5]. The radiative transfer equation for a plane-parallel atmosphere 
is given by 

cos 9 d -^ = -a(z)I„ + J u (z) (1) 

where (z) is the instantaneous radiant intensity that flows at each point in the medium per unit area, per 
unit of solid angle, at a given frequency a(z) is the extinction coefficient; and J is a source term. These 
last two quantities describe the loss/gain into the given direction. The angle 9 is the direction angle with 
respect to the vertical axis z with 9=0 when pointing upwards. 

In the general case, scattering into and from other directions can lead to both gains and losses to the 
intensity and are taken care by the terms a and J. For the microwave region, the scattering term is usually 
neglected [3]. If scattering is neglected, the only source term to consider is that due to local emission and the 
extinction coefficient reduces to the absorption coefficient a a . Assuming local thermodynamic equilibrium, 
each point can be characterized by a temperature T and from Kirchoff’s law we get 

J l ,(z)=<T a {z)B,A T'(z)) (2) 
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where B u (T) is the Planck function: 
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where h = 6.625 xlO -34 Js is Planck’s constant, c = 2.988x10 s m/s is the speed of light, and k =1.381x 10 : 
J/K. Equation (1) is a linear non-homogeneous first-order differential equation with solution 
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where I u (z 0 ) is the boundary condition, In the microwave region of the spectrum 

hv « kT 
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which results in (3) taking the form 
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(5) 


where A is the wavelength. This is known as the Raleigh-J eans approximation. From this expression, it is 
clear that in the microwave region the energy emitted is proportional to the physical temperature T. Another 
commonly used result from this relation is to define a scaling of the intensity 1„ as follows 



( 6 ) 


The quantity T b (v) is called the brightness temperature which is commonly used in the microwave retrieval 
literature instead of I n terms of brightness temperature and using (5), (4) takes the form 

T b {z) = T b (z 0 )e-M*")- 6 W + S e C (j J (7a ( 7 ) e -(ff(7)-«(*))sec^ 7 ^ 7 (7) 


where 

noo 

K z ) = J CaiOdt, 


( 8 ) 


is the optical thickness and co represents the top of the atmosphere (TOA). Here the dependency of all these 
quantities in frequency v is not shown for convenience 

For our purpose, it is of interest to solve this equation to obtain the brightness temperature that a satellite 
will measure at the top of the atmosphere when looking to the surface at an angle 0 off nadir. This will 
correspond toz 0 =0 surface(sfc) and z =oo in (7). The boundary term T b (0) is given by 


21(0) =e.T s -Kl - e)T d 


where T s is the surface temperature, T d is the downwelling radiation reflected by the surface back towards 
the satellite, and e is the surface emissivity. For the reflected component, it is assumed that the surface is 
a smooth, homogeneous, and isothermal so only the radiation in the specular direction 6 is accounted for. 
In our reference coordinates, the propagation angle for the downwelling radiation is n- 6. The downwelling 
radiation is obtained by solving the radiative transfer equation where integration is from the TOA to the 
surface with propagation angles - 0.This is taken care in (7) by setting z a =ooand z =0. The boundary 
term in this case is given by the cosmic microwave background emission with T c = 2.7 K. The resulting 
expression for the brightness temperature at the TOA at an angle 0 off nadir 

T b =[ eT s + (1 - e )T c e- 5 ( 0 ' se j= 8 e - 5 (°) sec ^ J T{z) [ - (1 - e)e -2(fi(0)-«(.)) sec J secO dz (9) 

Simulation of this expression if all quantities were known is a simple matter. However, the computation 
of the optical thickness and the absorption parameter requires the use of databases containing information 
about the spectral characteristics of atmospheric constituents such as HITRAN. I n our work, all optical 
depth computations were carried out using the FASE Radiative Transfer Code. 
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3 Atmospheric Remote Sensing 

The satellite instrument measures radiance that arrives into its field of view. The radiance that arrives is 
the sum of the radiance emitted and reflected by the surface, emitted and reflected by the atmosphere, and 
that scattered by the atmosphere into the field of view of the instrument. The relative contribution of each 
component depends on the region of the spectra seen by the instrument. Remote sensing of the surface takes 
advantage of those regions of the spectrum where the atmosphere is transparent or nearly so. I n the case of 
atmospheric remote sensing, the satellite sensor is looking at regions in the spectra where the atmosphere 
blocks the radiance emitted or reflected by the surface and therefore it receives that radiation that is emitted 
or reflected by the atmosphere. The interaction of electromagnetic waves with the atmosphere depends 
on the characteristics of the propagating wave (primarily its wavelength), the physical characteristics of 
the atmosphere and its constituents (pressure, temperature, density, absorbing gases, suspended particles). 
The mechanisms for interactions are: scattering, absorption, emission, and refraction. In regions where the 
atmospheric constituents characteristics are known or understood as in the 60 GHz oxygen and 183 GHz 
water vapor absorption lines, measurements of brightness temperature can be used to infer atmospheric 
properties of interest. The relation between atmospheric properties with brightness temperature is given by 
the radiative transfer equation (7). Therefore the problem of interest is to infer the atmospheric quantities 
of interest from measured brightness temperature by inversion of the radiative transfer equation. 


4 Temperature Retrieval Problem 

If the atmosphere strongly absorbs, most of the contribution to the measured brightness temperature will 
come from the atmosphere itself. In the case of the microwave region of the spectrum, scattering is negligible 
and the energy into the field of view of the sensor will come from atmospheric emission. Assuming that the 
satellite is looking at nadir (i.e. 0 =0), and that the surface temperature T s and emissivitye are known, (9) 
can be rewritten as 

pOO 

T b = T(z) I<(v,z)dz (lo) 

Jo 

where 

I<( v, z ) = [l -J- (l - e ) e -2(«(0MG»j e -«(») 

f h = T b - eT s + (1 - e) T c e" 5(0 j e~ 6 ^ 

If the absorber is uniformly mixed with a known concentration, as 0 2 , the quantity K(v, z) is known and the 
temperature profile T(z) could be retrieved by inverting (10). The function K(v,z) is called in the literature 
[3] the weighting function. From this point on, in our discussion we would not distinguish between T b and 
Tb in (10). 0 2 has several absorption lines around between 50 and 60 GHz. The Special Sensor Microwave 
Temperature-1 (SSM/T-1) sensor of the sensor suite of the DMSP satellite has 7 channels located in the 50 
to 60 GHz range used for temperature retrievals. A summary of the SSM/T-1 sensor characteristics is given 
in [6], In the 0 2 band, the shape of the weighting function is independent of the temperature making the 
inversion of (10) a linear inversion problem. 


4.1 Problem Discretization 

To numerically solve the temperature retrieval problem, first we discretize (10) by approximating the integral 
with a numerical integration formula. This results in the algebraic linear system of equation 


T fc = KT + e (U) 

where Ti6 7?. m and T 6 7 Z n are the brightness and atmospheric temperature vector; K 6 7 z mxn is the 
matrix of weighting functions; and e is an error term associated with measurement noise and the truncation 
error arising from the discretization of the integral equation. The number of measured brightness tempera- 
tures m, is usually smaller than the vertical resolution or number of temperature levels n to estimate. In the 
case of the SSM/T-1 sensor, there are m = 7 channels and normally n >20 temperature levels. Therefore 
the resulting algebraic linear system of equations (11) is under constrained (i.e. there are more unknowns 
than equations). 
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4.2 Regularization of the Discrete Problem 

The temperature retrieval problem is related to the solution of the linear system of equations (11). This 
problem has two major difficulties associated with it: (i) ill-conditioning due to the ill-posedness of the 
associated integral equation, and (ii) multiple solutions because of trying to estimate more temperature 
levels than measurements available. To overcome these difficulties, we will use the so called regularization 
methods. Regularization theory [1] transforms an ill-posed problem to a well-posed one, using a prion 
knowledge on the nature of the solution. Depending on the prior information, regularization techniques can 
be classified into two major groups: statistical and deterministic. 

4.2.1 Statistical Regularization 

In statistical regularization, prior statistical information is used to regularize the temperature retrieval 
problem. Our prior information in this case is the prior distribution of the temperature profiles p T (T) and 
the conditional distribution Pn/riTb)- According to Bayesian estimation theory [4], the best estimator T 
based on the brightness temperature observation T(, of the temperature profile T in a mean square sense is 
the conditional mean 

T=E(T/T 6 ) (12) 

We will refer to this estimator as the minimum mean square estimator (MMSE). The analytical determination 
of this function might be a very difficult task. In many instances, the estimator is constrained to be linear 
which results in the Linear Minimum Mean Squares Estimator (LMMSE)[4] 

T =T +A T ^ r , ; (T fc - T;,) (13) 

where T is the a priori mean of T, T ( , is the mean of T b given T, A T ,% is the cross covariance between 
T and T fc , and A TlJtTb is the conditional covariance of the brightness temperature. The LMMSE is easy to 
construct, since only the first and second order statistics are needed rather than their complete probability 
densities. Also, if T and T 6 are jointly Gaussian, the LMMSE is the optimal Bayesian MMSE. 

4.2.2 Deterministic Methods 

In this section, we will look at two regularization methods for ill-posed linear algebraic systems of equa- 
tions: Tikhonov regularization, and discrepancy principle regularization. Other methods are discussed in 
[2], Computation of the regularized solution was done using the MATLAB Regularization Toolbox presented 
in [2], 


Tikhonov’s Regularization 

One way to regularize (11) is computing T as the solution to the optimization problem 


T a 


: arg mm 
T£H"- 


||KT - Tb|| 2 + A 2 ||L (T — T 0 ) 


(14) 


where TO is a prior temperature profile estimate, || , \ is the 2-norm, and A is the regularization parameter. 

A key issue in this method is the selection of the regularization parameter A. The value used in the 
simulation results presented here was based on the L-curve method described in [2], The optimal value of A 
balances the prediction error :|KT- T 6 || with the regularization error ]|L(T- TO ) ||. 


Discrepancy Method 

Another possibility to regulate the temperature retrieval problem is by computing T as the solution to the 
quadratical ly constrained linear least squares problem 

T« = arg min ||L (T - T c )|| 2 (15) 

TeK” 

subject to ijKT - T b lj 2 <a 2 

where a plays the role of a regularization parameter. The solution to this problem can be made identical to 
Ta for a suitably chosen a [2]. We prefer to select a based on the measurement noise norm ||e||. 
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5 Simulation Experiments 

In this section, we present some simulation results that illustrate the use of deterministic methods to reg- 
ularize the temperature retrieval problem. For the simulation experiments, a 98 layer atmosphere based 
on the US Standard Atmosphere model was used. The optical thickness for each layer was computed us- 
ing theFASE code. The radiative transfer equation was numerically integrated using the trapezoidal rule. 
The implementation of Tikhonov’s and the discrepancy algorithm available in the Matlab™ Regularization 
Toolbox were used to compute the temperature retrievals. All computations with the exception of the optical 
thickness were done under the MATLAB™ environment. I n our simulations, the surface emissivit y e was 
set to 0.9 and the surface temperature T s set to 288.2 degrees Kelvin. Figures 1 and 2 show the results 
of applying Tikhonov’s regularization to the temperature retrieval problem with perfect measurements (no 
noise). The regularization parameter A was set 0.0044. The solid line is the retrieved temperature profile, 
the dashed line is the actual temperature profile, and the dash-dot line is the initial guess fed to the al- 
gorithm. We can see from this simulations that the algorithm was capable of estimating temperature up 
to 40 km. Beyond 40 km, the resulting estimate was identical to the initial guess. An important result is 
how the retrieval algorithm is capable of determining the height of tropopause corner when fed with initial 
guesses that have that corner at heights far from the actual height as shown in Figure 2. The importance of 
determining this peak comes from the fact that most important weather features are located at this region 
of the atmosphere. Also, the location of the tropopause peak serves as a figure of merit in evaluating the 
performance of temperature retrieval algorithms. We are not showing results for the algorithm based on the 
discrepancy principle regularization since they were similar to those of Tikhonov regularization.. 

Figures 3 and 4 show the performance of Tikhonov regularization under the presence of noisy data. The 
noise vector added to the measured brightness temperature has a normal distribution with zero mean and 
unit variance. The regularization parameter for this case was A =0,016. Notice that the noise causes the 
retrieved profile to be noisier with a maximum error in the first 20 km of 10 degrees Kelvin. Quite high 
compared to some retrieval methods that claim accuracies of 0.5 degrees Kelvin. However, the algorithm is 
still capable of determining the height of the tropopause. 

Figures 5 and 6 show the results for the retrievals in the noisy case computed using discrepancy principle 
regularization. The regularization parameter for this case was set at q = 3. That value corresponds to three 
standard deviations of the noise distribution. The resulting retrievals are smoother than those retrievals 
from Tikhonov regularization. The maximum error is in the neighborhood of 3 to 4 degrees Kelvin in the 
tropopause. The estimation of the location of the tropopause peak is also improved. 


6 Conclusions and Final Comments 

This paper presents some preliminary work in the application of regularization techniques to the I inear 
problem of atmospheric temperature retrievals from microwave radiometry. The proposed techniques were 
evaluated using simulated data. We used algorithms implemented in the MATLAB Regularization Toolbox 
developed by [2], The results obtained were quite encouraging. I n particular, being able to estimate the 
location of the tropopause corner of the temperature profile even with bad initial guess and noisy data is a 
result not previously observed with other algorithms. 
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Figure 1: Tikhonov Regularization with 10 km 
Tropopause. 



Figure 2: Tikhonov regularization with 20 km 
tropopause initial guess: noise free. 



Figure 3: Tikhonov regularization with 10 km 
tropopause initial guess: noisy case. 



Figure 4: Tikhonov regularization with 20 km 
tropopause initial guess: noisy case. 



Figure 5: Discrepancy regularization with 10 km 
tropopause initial guess: noisy case. 



Figure 6: Discrepancy regularization with 20 km 
tropopause initial guess: noisy case. 
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